# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# File Name          : run_paper_data_analysis.R
# Programmer Name    : Luis F. Campos
#                     luiscampos@g.harvard.edu
#
# Purpose            : This file contains code to run the analysis presented 
#					   in Section 7 of "Worth Weighting? How to Think 
#					   About and Use Sample Weights in Survey Experiments"
#
# Input              : 'survey_exp_toolkit.R', 
#					   'function_analysis_multilevel_treatment.R'
#					   'run_paper_data_analysis_sub.R'
#					   'Survey_Experiment_Data.Rda'
#
# Output             : final_results.csv
# 
# References         : None
#
#
# Platform           : R
# Version            : v3.3.0
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Replication Information:
#   We do expect differences when running this code on different systems 
#   (and with different seeds), but the differences all tend to be quite 
#   small and the general trends remain. If you would like to replicate 
#   the exact numbers and figures used in the article, in the readme you 
#   will find the session info. Please install and load all package 
#   versions (including the R version) to ensure the exact replication. 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #




set.seed(1)
# Set current working directory
# setwd('')

# load basic funcitons
source("./_survey_exp_toolkit.R")
# load analysis function
source("./_function_analysis.R")
# load experiments for analysis 
source("./_paper_analysis_labels.R")


out_dir = './'


library( plyr )
library( xtable )
library( reshape2 )
library(foreign, pos=4)
library(Matching)
library(car)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Number of Bootstrat samples for SE Estimates
#  Change to smaller number to speed things up
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

B = 1000

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Load Survey Experiment Data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #


load(file = 'Survey_Experiment_Data.Rda')


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	Virginia
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_Virginia = list()

for(i in 1:nrow(TMT_Virginia)){
	r = TMT_Virginia[i,]
	print(r)
	results_Virginia[[i]] = do.analysis(dat = dat_Virginia[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_Virginia, file = paste(out_dir, 'TMT_Virginia.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	CCES12_AGT_Matched_working
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_CCES12 = list()

for(i in 1:nrow(TMT_CCES12)){
	r = TMT_CCES12[i,]
	print(r)
	results_CCES12[[i]] = do.analysis(dat = dat_CCES12[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_CCES12, file = paste(out_dir, 'TMT_CCES12.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# YouGov
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_YouGov = list()

for(i in 1:nrow(TMT_YouGov)){
	r = TMT_YouGov[i,]
	print(r)
	results_YouGov[[i]] = do.analysis(dat = dat_YouGov[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_YouGov, file = paste(out_dir, 'TMT_YouGov.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 2010_CCES_UCB
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_CCES10_UCB = list()

for(i in 1:nrow(TMT_CCES10_UCB)){
	r = TMT_CCES10_UCB[i,]
	print(r)
	results_CCES10_UCB[[i]] = do.analysis(dat = dat_CCES10_UCB[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_CCES10_UCB, file = paste(out_dir, 'TMT_CCES10_UCB.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# CCES10
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_CCES10_EIT = list()

for(i in 1:nrow(TMT_CCES10_EIT)){
	r = TMT_CCES10_EIT[i,]
	print(r)
	results_CCES10_EIT[[i]] = do.analysis(dat = dat_CCES10_EIT[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_CCES10_EIT, file = paste(out_dir, 'TMT_CCES10_EIT.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	UCB_Follow
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_UCB_Follow = list()

for(i in 1:nrow(TMT_UCB_Follow)){
	r = TMT_UCB_Follow[i,]
	print(r)
	results_UCB_Follow[[i]] = do.analysis(dat = dat_UCB_Follow[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_UCB_Follow, file = paste(out_dir, 'TMT_UCB_Follow.Rda', sep = ''))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 2014_CCES_Matched
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
results_2014_CCES_matched = list()

for(i in 1:nrow(TMT_2014_CCES_matched)){
	r = TMT_2014_CCES_matched[i,]
	print(r)
	results_2014_CCES_matched[[i]] = do.analysis(dat = dat_2014_CCES_matched[[i]], w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = r[6], sub_level = r[7], B = B, return.all = TRUE)
}

save(results_2014_CCES_matched, file = paste(out_dir, 'TMT_2014_CCES_matched.Rda', sep = ''))




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	Load results and gather into one table
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
out_dir = './'

load(file = paste(out_dir, 'TMT_Virginia.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES12.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_YouGov.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES10_UCB.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_2014_CCES_matched.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES10_EIT.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_UCB_Follow.Rda', sep = ''))


results = c(results_Virginia, results_CCES12, results_YouGov, results_CCES10_UCB, results_2014_CCES_matched, results_CCES10_EIT, results_UCB_Follow)

nms = rep(c('Virginia', 'CCES12', 'YouGov', 'CCES10_UCB', 'CCES2014', 'CCES10_EIT', 'UCB_Follow'), c(length(results_Virginia), length(results_CCES12), length(results_YouGov), length(results_CCES10_UCB), length(results_2014_CCES_matched), length(results_CCES10_EIT), length(results_UCB_Follow)))


r.keep = c(1, 4, 7, 8)
morph = function(tab){
	nm = do.call('c', lapply(tab[r.keep,1], function(x) paste(x, c('.est', '.SE'), sep = '')))
	est = c(as.matrix(t(tab[r.keep,c(2, 3)])))
	names(est) = nm
	est
}

dats1 = lapply(results, function(x)x$dats)
res1 = lapply(results, function(x)morph(x$res))

res = data.frame(do.call('rbind', dats1), do.call('rbind', res1), survey = nms)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Neymanian Confidence Intervals and Significance
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
res$tau.sate.L = res[,"tau.sate.est"] - 1.96*res[,"tau.sate.SE"]
res$tau.sate.U = res[,"tau.sate.est"] + 1.96*res[,"tau.sate.SE"]
res$tau.hh.L = res[,"tau.hh.est"] - 1.96*res[,"tau.hh.SE"]
res$tau.hh.U = res[,"tau.hh.est"] + 1.96*res[,"tau.hh.SE"]
res$tau.ps.hh.L = res[,"tau.ps.hh.est"] - 1.96*res[,"tau.ps.hh.SE"]
res$tau.ps.hh.U = res[,"tau.ps.hh.est"] + 1.96*res[,"tau.ps.hh.SE"]

res$tau.sate.sig = ((res$tau.sate.L > 0) | (res$tau.sate.U < 0))
res$tau.hh.sig = ((res$tau.hh.L > 0) | (res$tau.hh.U < 0))
res$tau.ps.hh.sig = ((res$tau.ps.hh.L > 0) | (res$tau.ps.hh.U < 0))

res$delta = (res[,"tau.sate.est"]-res[,"tau.hh.est"])/res[,"tau.star1.SE"]



# This table contains all results from this analysis, it is called again and used 

write.csv(res, file = paste(out_dir, 'final_results.csv', sep = ''))

